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Abstract. We present a review on the mathematical methods used to theoretically 
study classical propagation and quantum transport in arrays of coupled photonic 
waveguides. We focus on analysing two types of binary photonic lattices where self- 
energies or couplings are alternated. For didactic reasons, we split the analysis in 
classical propagation and quantum transport but all methods can be implemented, 
mutatis mutandis, in any given case. On the classical side, we use coupled mode 
theory and present an operator approach to Floquet-Bloch theory in order to study 
the propagation of a classical electromagnetic field in two particular infinite binary 
lattices. On the quantum side, we study the transport of photons in equivalent finite 
and infinite binary lattices by couple mode theory and linear algebra methods involving 
orthogonal polynomials. Curiously the dynamics of finite size binary lattices can be 
expressed as roots and functions of Fibonacci polynomials. 
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1. Introduction 

The analogy between linear lattices and the atom-field interaction pQ or ion-laser 
interactions [2j [3] has been a fundamental step for the emulation, via classical 
interactions, of quantum mechanical systems. This is important due to both pure 
scientific interest and possible applications to quantum computing. In the latter, 
the properties of classical systems have been used to realize quantum computational 
operations by quantum-like systems and, in particular, it has been show how a 
controlled-NOT gate may be generated in nonhomogeneous optical fibers [I]. At the 
fundamental level, e.g., it has been possible to the emulate the most basic atom-field 
interaction, the Jaynes-Cummings model, theoretically [5] and experimentally [6] with 
arrays of photonic waveguides and, just to give another example, it has been proposed 
to model non-linear coherent states [7] in waveguide arrays [8]; linear coherent states 
have also been modelled in linear arrays of photonic waveguides [9]. 

In what follows we will take advantage of the simplest composite array of 
waveguides, i.e., binary arrays, to introduce the most basic mathematical methods used 
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to study photonic lattices. Section [2] gives an introduction to the two kinds of binary 
photonic lattices considered; those with either alternating self-energy or coupling. Then, 
we proceed to study the propagation of classical electromagnetic fields through binary 
lattices of infinite size by use of coupled mode theory for lattices with alternating self- 
energy and by an operator approach to Floquet-Bloch theory in the case of an array 
with alternating couplings. In section [4] we use coupled mode theory (with a twist 
given by the use of orthogonal polynomials) to solve a finite binary self-energies lattice 
in order to exemplify how the presented methods are valid, changing what needs to 
be changed, in all cases. In this section we also find the dynamics of a finite binary 
couplings array by basic matrix theory methods. Curiously, we find that the dynamics 
for both finite size binary lattices can be written in terms of Fibonacci polynomials 
evaluated at the roots of a Fibonacci polynomial which order is related to the size of 
the system. In section [5] we introduce some quantities of interest when studying photon 
transport in arrays of photonic waveguides. We show that initial states with a Gaussian 
distribution of amplitudes and linear coherent states, that is, initial states with Poisson- 
like distributions, reconstruct in binary lattices. Finally, we present our conclusions. 

2. Binary photonic lattices 




Figure 1. (Color online) (a) Photonic binary super-lattice where identical waveguides 
alternate different nearest neighbor coupling, (b) Photonic binary super-lattice where 
homogeneously coupled waveguides alternate diffraction index. 

A binary photonic super-lattice is composed by the repetition of a primitive unit 
cell where two different elements are characterized by one parameter. One type of such 
binary waveguide arrays, shown in figure [IJa), is composed by identical waveguides 
where nearest neighbour coupling between them alternates between two values; hereby, 
we will refer to this type as binary coupling (BC) lattice. In the other type, shown 
in figure [l](b), waveguides with two different refraction indexes alternate position and 
are homogeneously coupled; which we will call binary index (BI) lattice from here on. 
The propagation of a classical light field in an infinite BC lattice is ruled by the scalar 
differential equation set for the field amplitude at the jth waveguide, 

- id z Ej = nEj + g (E j+1 + + (-l) j S (E j+1 - E^) , j = -oo, . . . , oo, (1) 

where the constant n is the refraction index of each waveguide and a total coupling has 
been defined as go = (gi + #2) with the difference as 5 = (g\ — #2) /2. We have used the 
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shorthand notation d x to express the partial derivative with respect to x. An infinite 
BI lattice it is given by, 

- id g Ej = [n + (-l) j e\ E 3 + g {E j+1 + Ej- X ) , J = -oo, . . . , oo, (2) 

where a base refraction index has been defined as n Q = (ni + n 2 )/2 such that it is 
halfway between the refraction index at each waveguide, n\ and n 2 , i.e., e = \n\ — 712I/2, 
and the tight binding coupling is given by the real constant g. Infinite [TOj [III 121 [13] 
and semi-infinite [2] BC and BI lattices, with the addition of a non-linearity to each 
waveguide, are well known in the non-linear optics community where the existence and 
stability of diverse types of solitons have been studied along the years. 

In linear optics, it is possible to borrow from condensed matter theory a Floquet- 
Bloch result for quasiparticle motion on a chain in order to solve the evolution of a 
classical field, or a single photon, through these lattices [T5]; e.g., the dispersion relation 
for a BC lattice, shown in figure [2^ a), 

^ = 4^ 2 + (^-^)cos 2 0], (3) 

results of the infinite BC lattice has been used to discuss the existence of two propagation 
modes with opposite transverse velocities [TO]- For a BI lattice the dispersion relation 
is, 

fiJ = /3 2 + 4cos 2 0, P = e/9 (4) 

The band gap structure of this dispersion relation, shown in figure |2^b), near the edge 
of the Brillouin zone, 4>b = 7r/2, suggests the use of Bloch waves close to 0# to simulate 
one-dimensional Dirac equations [TTJ . Under such a condition, a photonic analogue of 
zitterbewegung has been theoretically proposed [18J and experimentally realized 
more complex approaches of such a classical simulation include atoms in bichromatic 
optical lattices pTJ| [21] . 
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Figure 2. (Color online) Dispersion relation for infinite (a) BC with go — 35 and (b) 
BI lattices. 

3. Classical electromagnetic field propagation 

3.1. Coupled mode theory 

We are going to implement a coupled mode theory method on the infinite BC lattice 
described by (flj. For starters, let us define the field amplitudes with a rotation 
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proportional to the refractive index function, n , i.e., Ej — > e m ° z Ej] in this way, we 
can get rid of the term involving the refractive index function, 

- \d z Ej = g (E j+1 + Ej^) + (-l) j S (E j+1 - E^) , j = -oo, . . . , oo. (5) 

Now, we can define two auxiliary proper modes, say A k = e lQqZ E 2 k and B k = e lQqZ E 2 k-i, 
such that our differential set is now given by a coupled modes set 

(go + 5)A k + (g - 5)A k+1 = n q B k , (6) 
(go + S)B k + (g - = Sl q A k . (7) 

It is trivial to find a three term recurrence relation for the mode B from this coupled 
set of equations, 

A k ^ l + A k+l nl-2(gl + 5 2 ) 



A k g 2 - 5< 



(8) 



solved by 



A k = e 1 ^, (9) 

2 e i(2fc+l)0 

B k = — - [g cos + i<5 sin 0, ] (10) 



2 cos 20= * ™ (11) 



under the restriction 

% - (4o 2 + 5 2 

From the expression above, we straightforwardly recover the dispersion relation 

^ = 4[5 2 +(^-5 2 )cos 2 0]. (12) 
Notice that we can also write the normal mode components and dispersion relation as, 

j even, 
j odd, 



E ( ? ] = / ' (13) 



J I I go cos <j>+i8 sin 

g cos <j>— i<5 sin </> 

= \gcos<f) + i<5sin0| 2 . (14) 

Under this treatment we can recover the field amplitude propagation at the jth 
waveguide for a field entering the lattice at the rath waveguide, 

E M 1 r d(f> e i( m - 3 >5W ) (15) 
2tt 

1 m — j even, 



( ?r^r/ )" j even, m odd, (16) 
f goccej+iBffof Y j odd, m even, 

k V 9 cos <p— id sm (p I 

and we have found all the necessary information to study the propagation of any given 
initial classical field impinging a BC lattice. 
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3.2. An operator approach to Floquet-Bloch theory 



It is possible to use an operator approach to find Floquet-Bloch waves. In order to 
demonstrate it, let us write the differential equation set for an infinite BI lattice ^ in 
matrix form as 



- \d z E = HE, 
H=(-lff3 + V + V\ 



with solution 



E(z) 



£7(0), 



E(z) = J2 E MEj 



(17) 
(18) 

(19) 



J=-oo 



The expression (17) is identical to the whole differential set Q up to a unitary 
transformation equivalent to a change into a frame rotating at a frequency proportional 
to the refractive index function no and a change of units in terms of the homogeneous 
coupling g. We used the unitary ladder operators, VW = VV^ = 1, defined as 



V\k) = \k-1), V^\k) = \k + 1), h\k) = k\k). 



(20) 



Where we substituted E^ — > \k), i.e., the state \k) represents the field localized at the 
fcth waveguide. These ladder operators fulfil the commutation relations 



n,V 



-V, 



n 



0. 



By using this operator representation, we can define a phase state as 

oo 

10) = eifc ^>> 



(21) 



(22) 



which is the Fourier transform of the localized field. In other words, these operators 
allow us to do Floquet-Bloch theory, 

' 1 d0e-^|0). (23) 



2tt 



In this phase state basis, the action of the operator Hbi over some useful phase states 
reduces to: 



H 2 



+ 71") + 2 cos 



(/3 2 + 4 



cos 



H 2 \(f> + ir) = (/3 2 + 4 cos 2 



+ 7T>, 



(24) 
(25) 
(26) 



From the last two equations it is possible to infer the dispersion relation, 

ft J = /3 2 + 4cos 2 0. (27) 

In order to recover the rest of the information given by Floquet-Bloch theory, let us 
calculate the evolution of our phase state, 



AHz\ 



E 



\2k 



{2k)\ 



+ 



,2fc+l 



(2k + l)\ 



H 2 ) H 



cos fl^z] 



+ i 



sin Vt^z 



{P\4> + tt) +2 cos <P\<P)). 



(28) 
(29) 
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Now, if the field started at the mth waveguide, |^(0)) = |m), one can write the field 
amplitude at the jih waveguide as Ej = (j\m(t)} and, by use of (m\4>) = e im< ^, recover 
the field amplitude evolution from standard Floquet-Bloch theory 

E (m) 1 r d(j) e i( i - m) ^(m) ) (30) 



— 7T 



Mm) n , . ^CQS^ + ^-i)^ 

£• = cos i 2^2; + i — smil^z. (31) 

V sfy J 

4. Photon transport 

In some cases one has to consider the evolution of a quantum field in an array of 
coupled photonic waveguides; in a general case, it could be arrays of coupled microring 
resonators, coupled cavities, or capacitively coupled strip-line resonators. When such a 
system is presented it is usually described by a Hamiltonian which, in the case of finite 
BC and BI lattices of size N, is given by 

Af-1 

Hbc =^2 [oo- {-1Y8] + > ( 32 ) 

3=0 

N N-1 

H /,. 

3=0 j=0 

where we have moved into an adequate reference frame rotating at a frequency defined 
by the refractive index function no and set units in terms of h and hg, respectively. The 
operator a\ (a&) creates (annihilates) a photon in the kth. waveguide. 

In the Heisenberg picture, the equations of motion for these systems are 

BC : id t CLj = g (a j+l + Oj_i) + (-1) J (a j+1 - , (34) 
BI : idtdj = (— 1) J (3aj + hj-i + a,+i, (35) 

these sets with j = 0, 1, . . . , N, are the finite equivalent, by making the substitutions 
fij — > —Ej alongside t — > z , of the differential sets ruling the propagation of a classical 
field through the corresponding binary photonic lattice in ^ and (17), respectively. 



JV JV — 1 



4-1. Coupled mode theory 

Let us start with the finite BC lattice and introduce the transformation 

f _ e iffo Ef^fatf+i+zfe+i) ; (36) 

such that a general states is defined by = T\<f>) and leads to the Schrrodinger equation 
in units of hS 

id t \4>) = H\4>) 1 (37) 

Af-1 

H = -(-^ (Vl+i + a J a i+i) > ( 38 ) 

j=0 
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with equations of motion in the Heisenberg picture 

idtdj = {-I) 3 (a j+1 - a^i) , j = 0, 1, . . . , N. (39) 

Now, let us define the coupled modes as a 2 fc(t) = — ie n *&2fc+i and a 2 k+i(t) = e nt b 2 k+2- 
Notice that we have used Qt in the argument of the exponential instead of iflt, thus we 
are looking for a purely imaginary Q. This allows us to write an eigenequation in the 
form, 

N-l 

Mb = Qb, M id = 8ij- x + 6i- ltj , b =J2 c A+i, ( 4 °) 

3=0 

where it is possible to recover a recurrence relation for the coefficients of the normal 
modes, 

c 2 = flci, (41) 
c 3 = Qc 2 + ci, (42) 

; (43) 

Qc N + Qc N ^ = 0, (44) 

that is solved by Fibonacci polynomials, 

Ci = ^(n). (45) 

The last of the recurrence relations gives a boundary condition that translates into the 
expression Fn+i(Q) = that is solved by [22J 

m = l ±2isin f^ eveniV ' (46) 

V ; 1 i^sin^ odd N, V ; 

and we recover the purely imaginary eigenvalue that we were looking for. So, the normal 
modes dk = Ylj=o iF2j+i(i\k)a>2j + F2j+2(^k)^2j+i (up to a normalization constant) 
with eigenvalues = —ifl(k), diagonalize the Hamiltonian H = J^f^d ^Ac\. and we 
have found the dynamics of the system. It is trivial to go back into the original frame. 

4-2. A linear algebra approach 

Now, let us consider a finite BC lattice. Heisenberg equations of motion for this system 
may be written as the matrix differential equation, 

d t a = -iMa, Ma = (-1)^ + <^-i + S i>j+1 , a(t) = e" iA/ 'a(0), (47) 

The solution to this system may be found in a number of ways [23] and we do it by 
finding the system, {V, A}, of the matrix M = VAV* 1 [24J; where the eigenvalues matrix 
A is a matrix which diagonal elements are the N eigenvalues of the matrix M, {Xj}, 
and the eigenvector matrix V is a matrix which rows, Vj, are the N eigenvectors of M. 

The characteristic polynomial, p^, of tridiagonal matrix M is found via the method 
of minors (25 1 as 



(-1)^F N+1 (^P^V) iVeven, 

PnW = ) (-1) (7V " 1)/2 ^^^(v^ 3 ^) N odd, ^ 
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(-l)-N/ 2 b N / 2 ((3 2 - A 2 ) iVeven, 
(_l)(iv-i)/2 (/ 3_ A) B is - 1)/2 (p - \*) iVodd, l4yj 

where F n (x) is the nth Fibonacci polynomial [20], ^n(^) and B n (x) are nth Morgan- 
Voyce polynomials |27j . The roots of Fn+i{x) = are well known [22] and yield the 
eigenvalues, 



= - J^+W^ ) , < AT/2, 

^ +4c0S 2(j£_) 3 '>JV/2. 

These proper values, alongside (M — Ajl) = 0, deliver recurrence relations fulfilled by 
eigenvector components, 

V jlk = r=== , j,k = 0,...,N-l, (51) 

where orthogonal polynomials, defined as: 

f (-l)^ Ffc+1 (2i|cos^|) fceven, 
1 -i(-l) (fc+1)/2 ^£|^^ +1 (2i|cos^|) k odd, ( 52 ) 



(-l) fc/2 & fc/2 (-4cos 2 ^) fceven, 
(_l)(fe + i)/2 (/3 _ A .) S(fc _ 1)/a (_ 4cos 2 _^_) fc odd . 



(53) 



Thus the Hamiltonian is diagonalized, H = ^2f = Q XjCjCj, in terms of delocalized normal 
modes c = V^a, giving a time evolution a(t) = e~ lMt a(0) = Ve~ lAt V"^a(0). 



5. Quantities of interest 

For the sake of space and simplicity we will only discuss typical quantities of interest 
when studying the propagation of a quantum electromagnetic field in an array of 
photonic waveguide lattices. It is trivial to show that the time evolution of an arbitrary 
initial state can be obtained from the evolution of the annihilation operator found in 
the last section. In order to generalize, let us write such an evolution as 

N-l 

% = S^W^( )> (54) 

A:=0 

where it is trivial to construct the matrix U(t) from the methods given in the last 
section. 

The most basic and visually appealing quantity of interest is the transport of m 
photons impinging a single- waveguide at time zero, 

|^(0)) = _Lat»»(o)|o), p = 0, . . . , N - 1. (55) 

The mean photon number at the gth waveguide for such an initial state after propagation 
is given by 

(n q ) p = (^(0)|a(*)ta(i) ? |^(0)> = m\U M (t)\ 2 . (56) 
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In the case of single photon input, m = 1, the expression for the mean photon number 
is identical to the normalized intensity at the qih waveguide for a propagating classical 
electromagnetic field. 

Superposition states of single photons are also interesting, 

N-l N-l 

|^(o)) = 5>^(o)|o>, Yl N 2 = x > ( 5T ) 

3=0 3=0 

and we can easily calculate the probabilities of finding the photon in the gth waveguide, 



N-l 



2 



(5f 



a 3 U P,l( t ) ■ 
3=0 

Here, we are interested in two peculiar single-photon distributions. One that we will 
call a Gaussian-like input distribution, 

JV-l 

|^,(0)) = e- k V {2 < )+iqk i /2 a)\Q), kj =j- \N/2] , (59) 

3=0 

which could be thought as the result of a Gaussian beam whose intensity peak is aligned 
with the center of the lattice impinging parallel, q = 0, or at an angle, q ^ 0, with the 
propagation axis [TSJ [19]. The other is similar to a quantum coherent state and is 
produced by a single-photon entering a Glauber-Fock lattice at the zeroth waveguide 
[SEIGI], we will call this a Poisson-like input distribution, 

N-l j 

| a ) = e -M"/>£ « & t|o>. (60) 



3=0 yJ 

In these Gaussian- and Poisson-like distributions, the maximum probability to find the 
single-photon is given at the |~iV/2]th and |~|ct| 2 ]th waveguides (where \x\ has been used 
to express x rounded up), in that order, and the variance is given by w and \N/2], 
respectively. Both initial states have a momentum in the direction perpendicular to 
propagation given by q and Im(a) in the sense that the center of mass for an input with 
q 7^ or Im(a) 7^ will drift to the right or left side of the lattice depending on the 
sign of q or Im(a). Notice that a Gaussian-like distribution will always be symmetric 
with respect to the center of the lattice, which is not the case for the coherent-like 
distribution. Propagation of Gaussian-like input distributions have been studied in 
infinite BC and BI lattices. In the former, the lattice splits the q = input in two 
components propagating in opposite directions [16]. In the latter, BI lattices allow the 
clasical simulation of the Dirac equation when the parameter q is close to the edge of 
the Brillouin zone [HJ [19] . To our knowledge there exist no record in the literature of 
propagation of Poisson-like input distributions through finite binary lattices. Figure [3] 
shows the mean photon probability for two kinds of superposition states impinging a BI 
and BC lattice, respectively. It is possible to see that the states reconstruct periodically 
in both cases through the fidelity function 

F(t) = \(Ma)\Mt))\ 2 - (6i) 
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Figure 3. (Color online). Time evolution of (a,c,e,f) the mean photon number 
for (a,e) Gaussian- and (c,g) Poisson-like input distributions with parameters 
{wo = 7, q = 0.557r} and a — \/50, in that order, and (b,d,f,h) their respective fidelities 
in a (a-d) BC and (e-h) BI with (3 = 0.5 lattice of size TV = 101. Time in units of 5 for 
BC and g for BI lattice. 



It is also possible to study the time evolution of the center of mass, 

AT-1 

jcm = 22 k(%a>k), (62) 

k=0 

for Gausian- and Poisson-like input distributions with complex parameters. Time 
evolution of the Fidelity for Gaussian-like distributions heavily fluctuates with values 
below 1/2, see top-right insets in figure |4^b,f), and its center of mass localizes at the 
central waveguide after a long time; after t = 7500(7 i n the case shown in figure |4](f). 
Only Poisson-like distributions partially reconstruct, see top-right insets in figure |4][d,h.). 
Of course, well-known results regarding classical simulation of Dirac equation |T8j are 
reproduced; for instance optical zitterbewegun by using a Gaussian-like distribution, 
figure [4j^f ) . Notice that the center of mass of the coherent-like distribution, figure Qh), 
wobbles with an almost negligible amplitude compared to that of the Gaussian-like 
distribution, figure PJf), and it appears to reconstruct with low fidelity after the second 
edge reflection; top-right inset in figure |l](h) . 

Another interesting set of states are product states, 

k 

1^,(0)) =n«l|0), (63) 
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Figure 4. (Color online). Time evolution of (a,c,e,f) the mean photon number for 
(a,e) Gaussian- and (c,g) Poisson-like input distributions with complex parameters 
{ojq = 7, q = 0.557t} and a — V50 — 0.557T + i\/0.55ir, in that order, and (b,d,f,h) 
their respective center of mass (with top-right insets showing the fidelities and top- 
bottom a longer-time evolution of the center of mass) in a (a-d) BC and (e-h) BI with 
/3 = 0.5 lattice of size TV = 101. Time in units of S for BC and g for BI lattice. 



where x — (xi, . . . , x k ) with X{ e [0, N — 1] and Xi ^ Xj for any i ^ j. The evolution of 
the photon number at the gth waveguide for product states is given by 

k 
3=1 

An example of these states is the two-photon product state given by 

|^ 8 (0))=a]al|0), (65) 
that gives a mean photon number evolution and two-photon correlation function 

(hq)ps = \Uj, q \ 2 + \U kt g\ 2 , (66) 

^ } = \U P ,jU^ k + U Pyk U qj \ 2 . (67) 

For the sake of space we will finish this section presenting how to deal with NOON 
states. A higher order NOON state, their mean-photon evolution at the gth waveguide, 
and two-photon correlation are given by the expression, 

\m) = -U U m + ^1™) |0>, m = 2, 3, . . . , (68) 
9.\/m,' V / 
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777/ 

(n q ) = j (P^ + \U k , q \ 2 ) , (69) 
T P , q = \U pJ U qJ \ 2 + \U p , k U q>k \ 2 + 2Re (e^U^U^U^U^) . (70) 

6. Conclusion 

We presented a review on mathematical methods used to study infinite and finite 
photonic lattices and the quantities of interest for quantized electromagnetic fields. Each 
method is presented in a particular context but any of the methods can be used in all case 
by making the necessary alterations. Alongside this review, we introduced a novel result, 
up to our knowledge, in the form of the exact spectra and proper modes for BC and 
BI lattices of size N (BI latices include the single-type lattice when the characteristic 
parameter is zero) in terms of the roots of the N + 1 Fibonacci polynomial and the 
Fibonacci polynomials evaluated at these roots, in that order. In order to illustrate our 
results we chose to focus on the evolution of multiple-waveguide single-photon inputs 
in the form of what we called Gaussian- and Poisson-like distributions impinging odd 
lattices; these distributions are related to Gaussian beams and to the output from 
Glauber-Fock lattices, respectively. Due to their spectral decompositions, Gaussian- and 
Poisson-like distributions with real parameters partially reconstruct quasi-periodically 
when they impinge a binary super-lattice with their intensity peak aligned with the 
middle of the lattice. 
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